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Using the model of Nambu and Jona-Lasinio to provide a microscopic description of both the 
structure of the nucleon and of the quark to hadron elementary fragmentation functions, we inves- 
tigate the transverse momentum dependence of the unpolarized quark distributions in the nucleon 
and of the quark to pion and kaon fragmentation functions. The transverse momentum dependence 
of the fragmentation functions is determined within a Monte Carlo framework, with the notable 
result that the average P\ of the produced kaons is significantly larger than that of the pions. We 
also find that (P±) has a sizable z dependence, in contrast with the naive Gaussian ansatz for the 
fragmentation functions. Diquark correlations in the nucleon give rise to a nontrivial flavor depen- 
dence in the unpolarized transverse-momentum-dependent quark distribution functions. The (fey) 
of the quarks in the nucleon are also found to have a sizable x dependence. Finally, these results are 
used as input to a Monte Carlo event generator for semi-inclusive deep inelastic scattering (SIDIS), 
which is used to determine the average transverse momentum squared of the produced hadrons 
measured in SIDIS, namely, (P T ). Again, we find that the average P T of the produced kaons in 
SIDIS is significantly larger than that of the pions and in each case {Pt} has a sizable z dependence. 

PACS numbers: 13.60.Hb, 13.60.Le, 13.87.Fh, 12.39.Ki 
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I. INTRODUCTION 

Semi-inclusive deep inelastic scattering (SIDIS) has a 
very rich structure which provides a wealth of observ- 
ables far in excess of the familiar inclusive deep inelastic 
scattering (DIS). The 2-dimensional picture of a target 
provided by SIDIS promises many new insights into nu- 
cleon and nuclear structure [DS]- F° r example, it has 
been realized that SIDIS may shed light on the angular 
momentum structure of the proton in terms of the spin 
and orbital angular momentum of its quarks and glu- 
ons [S]-[7]. It will also provide new information on the 
in-medium modification of bound nucleons and deepen 
our understanding of QCD itself PHI]- The study of the 
transverse momentum distribution of hadrons produced 
in SIDIS [IHH I5HTT] is characterized by determining the 
transverse-momentum-dependent (TMD) parton distri- 
bution functions (PDFs) and the TMD fragmentation 
functions. 

Early theoretical models of the fragmentation func- 
tions have been constructed in Refs. [T^l - fT!)] and more 
recently the development of the NJL-jet model [17] has 
provided a framework which automatically satisfies the 
relevant sum rules. Lattice QCD studies of TMD PDFs 
are presented in Ref. [18] and the QCD evolution of TMD 
PDFs is discussed in Ref. [19] . Extensive phenomenologi- 
cal data analysis of transverse momentum in distribution 



and fragmentation processes was presented in Ref. |20) . 
Considerable experimental work has already been car- 
ried out at JLab [2TH25], HERMES [26H2S] and COM- 
PASS [2"9"rl3"T] . while for an overview of the future per- 
spectives for this field we refer to the recent review by 
Anselmino et al. [52"] . 

In this work, we present the first microscopic calcu- 
lation of the spin-independent TMD quark distribution 
functions in the nucleon and the TMD quark to pion 
and kaon fragmentation functions, where none of the pa- 
rameters are adjusted to TMD data. The underlying 
theoretical framework is the Nambu- Jona-Lasinio (NJL) 
model [33J OH] . While this certainly represents a simpli- 
fication of QCD, it has many desirable properties. For 
example, it is covariant and respects the chiral symmetry 
of QCD, including its dynamical breaking. Moreover, it 
describes the spin and flavor dependence of the nucleon 
PDFs, as well as their modification in-medium [35H37] . It 
also produces transversity quark distributions 38 which 
are in good agreement with the empirical distributions 
extracted by Anselmino et al. |39| . 

For the present purpose the recent developments in the 
NJL-jet model [TTl [40], [41] , which provides a quark-jet de- 
scription of the fragmentation process using elementary 
fragmentation functions calculated within the standard 
NJL model, are also critical. This framework provides a 
good description of the parametrizations of experimen- 
tal data and has been extended to include vector me- 
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FIG. 1. Illustration of the three-dimensional kinematics of 
SIDIS. The photon momentum defines the z axis and the 
struck quark has initial transverse momentum fcx in the nu- 
cleon, with respect to the z axis. 



son resonances and nucleons as fragmentation channels. 
The use of Monte Carlo methods to calculate these frag- 
mentation functions has also been implemented and that 
development allows us to address a wider array of pro- 
cesses within the model, including physical cross-section 
calculations. 

In Sec. |TlJ we present the general formalism for de- 
scribing transverse momentum distributions in SIDIS, in- 
cluding the quark-jet model originally proposed by Field 
and Feynman. The calculation of the elementary, unintc- 
grated fragmentation functions in the NJL model, which 
are the input to the jet model which describes the TMD 
fragmentation functions in quark hadronization, is ex- 
plained in Sec. 
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Our model for the TMD quark dis- 
tribution functions in the nucleon is outlined in Sec. IV 
where we also present results for the TMD PDFs. Re- 
sults for the TMD fragmentation functions are discussed 
in Sec. [V] and the average transverse momentum in the 
SIDIS process, determined using our Monte Carlo event 
generator, is discussed in Sec. |VI| Finally, Sec. |VII| con- 
tains a summary and outlook. 



II. TRANSVERSE MOMENTUM IN THE 
NJL-JET MODEL 

The kinematics of semi-inclusive hadron production, 
IN — > I'hX, is illustrated schematically in Fig. [TJ where 
a lepton with momentum I scatters on a target, by emit- 
ting a virtual photon with momentum q that hits a quark 
with initial momentum k. As usual, the z axis is chosen 
to coincide with the direction of the photon's momentum, 
where the target has its momentum in the negative z di- 
rection. The transverse momenta in the process - labeled 
with a subscript T - are defined with respect to this z 
axis, so that the photon and target have no transverse 
momentum component (jN collinear kinematics). The 
angle between the lepton scattering plane and the quark 




FIG. 2. Illustration of the kinematics of SIDIS, where the 
final transverse momentum of the produced hadron with re- 
spect to the z axis is denoted by Pt, which is related to the 
initial quark transverse momentum in the nucleon fcj and that 
generated in the fragmentation process P± by Eq. jl]). 



scattering plane is denoted as ip. We allow for the struck 
quark in the target to carry a transverse momentum fey. 
Some of this transverse momentum is then transferred to 
the hadrons emitted by the quark. 

The kinematics of the quark fragmentation process is 
depicted in Fig. [2j The emitted hadron h carries a trans- 
verse momentum Pt with respect to the z axis which 
can be decomposed into two contributions. First, the 
quark transfers a fraction of its transverse momentum 
kx to the hadron and second the hadron also acquires 
a momentum transverse to the direction of the quark's 
momentum, P±. Up to corrections of order 0(k T /Q 2 ), 
the following relation holds [32]: 



T = *~±_ 



kn 



(1) 



This relation allows one to probe the quark transverse 
momentum inside a nucleon by measuring the z de- 
pendence of the emitted hadron's transverse momen- 
tum (P T ), provided (Pj_) is independent of z. How- 
ever, in the NJL-jet model framework we find that (P]_) 
is strongly z dependent and this z dependence is also 
observed at COMPASS [3T]. A recent analysis of the 
HERMES data [27] was performed in Ref. [20], where 
a Gaussian ansatz for the TMD quark distribution and 
fragmentation functions was assumed and an average was 
performed over the quark flavor and type of hadron de- 
tected. Using a fit region of 0.2 < z < 0.7, they extracted 
the following results for the average transverse momen- 
tum squared [20 : 



(k T ) = 0.38 ± 0.06 GeV 2 
(Pi) = 0.16 ±0.01 GeV 2 



(2) 
(3) 



The latest iteration of the NJL-jet model [UJ em- 
ploys Monte Carlo simulations to calculate the integrated 
quark fragmentation functions. It assumes that the ini- 
tial high energy quark emits hadrons in a cascade-like 
process, schematically depicted in Fig. [3] At every emis- 
sion vertex we choose the type of emitted hadron h and 
its fraction of the light-cone momentum z of the frag- 
menting quark, by randomly sampling the correspond- 
ing elementary quark fragmentation (splitting) functions, 



FIG. 3. NJL-jet model including transverse momentum. 

d,q(z), that are calculated within the NJL model. In each 
elementary fragmentation process we record the flavor of 
the initial and final quarks and the type of the emitted 
hadron, we also note the light-cone momentum fraction of 
the initial quark transferred to the hadron and that left to 
the final quark. The fragmentation chain is stopped after 
the quark has emitted a predefined number of hadrons, 
N Links- We repeat the calculation Nsi ms times, with 
the same initial quark flavor, q, until we have sufficient 
statistics for the emitted hadrons. The fragmentation 
functions are then extracted by calculating the average 
number of hadrons of type h, with light-cone momentum 
fraction z to z+Az, which we denote by (N q (z, z + Az)). 
The fragmentation function in the domain [z, z + Az] is 
then given by 

v J2n N^(z,z + Az) 
D h q {z)Az = (N h q (z, z + Az)) = qK '. 

■L* Sims 

(4) 

In this work, we extend the NJL-jet model to include 
the transverse momentum dependence of the emitted 
hadrons in the fragmentation process. This is achieved 
by using TMD elementary quark fragmentation functions 
at the hadron emission vertices and by keeping track of 
the transverse momenta of all the particles in the process. 
Our goal is to calculate the TMD fragmentation function, 
D q (z, Pj_), using its probabilistic interpretation. That is, 
the probability of a quark q to emit a hadron h with a 
fraction z of its light-cone momentum and a transverse 
momentum P± is given by D q (z, P]_) dz d 2 P±. 

We calculate elementary (one-step) TMD splitting 
functions, W^{z,p\), using the NJL model, where p± de- 
notes the transverse component of the hadron's momen- 
tum with respect to the parent quark, as illustrated in 
Figs. [3] and [4| In each step of the Monte Carlo simulation 
of the quark cascade emission, we randomly sample the 
type, the light-cone momentum fraction, z, and the trans- 
verse momentum, p±, of the emitted hadron using as the 
probability distribution the elementary TMD splitting 
functions of the quark, where the elementary probabil- 
ity is d(z,p 2 _^) dz d 2 p±. Schematically, the quark emission 
process is depicted in Fig. |4j where the z' axis denotes 
the direction of the original parent quark's 3-momcntum. 
The vectors k and k' denote the 3-momentum of an arbi- 
trary quark in the cascade chain before and after hadron 
emission with transverse components k± and k' ± , respec- 



FIG. 4. Quark elementary fragmentation kinematics, for an 
arbitrary hadron emission in the cascade chain. The z' axis 
is defined by the direction of the 3-momentum of the original 
parent quark. 

tively. The emitted hadron's momentum is labeled by Ph, 
where its transverse component with respect to k and the 
z' axis is denoted by p± and P± , respectively. P± is ob- 
tained using the relation P± = p± + z k± , analogous to 
that in Eq. (fTl) . The recoil transverse momentum of the 
final quark, is calculated from momentum conserva- 
tion in the transverse plane, namely, 

k ± = P ± + k' ± . (5) 

The TMD fragmentation function is then calculated after 
the trivial integration over the polar angle of P± in the 
transverse plane, that is 

D h q (z, Pi) AzttAP'I = (N%{z, z + Az, P\.Pl + AP'l) 

^N 3ims N£(z, z + Az, Pi, Pi + API) 
= N ■ (6) 

Jv Sims 

The model can easily accommodate the initial trans- 
verse momentum of the quark, for example, with re- 
spect to the direction of the virtual photon in SIDIS (see 
Fig.[T]). Our goal is to describe the average transverse mo- 
mentum of the hadrons produced in different reactions. 
The differential cross section for SIDIS up to terms of 
order 0(k 2 -./Q 2 ) can be written as [32] 

1 q 

^J24D h q (x,z,P 2 ), (7) 
1 

where Pj_ and Pt are related by Eq. ([T]) and q(x, k T ) 
are the TMD quark distribution functions of the target. 
Thus for SIDIS, we can use the TMD quark distribution 
functions to randomly sample the initial transverse mo- 
mentum of the quark to calculate the relevant number 
density of the produced hadrons. In this work, we use 
TMD valence quark distributions in the nucleon - calcu- 
lated within the NJL model - to determine the average 
transverse momentum of the produced hadrons with re- 
spect to the direction of the virtual photon, that is (P T ), 
by calculating the corresponding probability densities 
D q (x,z,P T ), using an expression analogous to Eq. d6j). 
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In this way, we obtain a self-consistent description of the 
entire process in the regime where the virtual photon 
samples the valence quark component of the target, that 
is, when the struck quark has x > 0.3. The type of target 
and the allowed range of x in the Monte Carlo simula- 
tion can be matched to those measured in any particular 
experiment. 

In this article, we only consider the production of pseu- 
doscalar mesons, that is, the pions and kaons, as a first 
step in determining the TMD fragmentation functions. 
Eventually, we will also include the vector mesons and 
nucleon-antinucleon channels, as done for the integrated 
fragmentation functions in Ref. [4"T] . 

III. ELEMENTARY TMD FRAGMENTATION 
FUNCTIONS 

In this section, we evaluate the "elementary" fragmen- 
tation functions of quarks to hadrons as a "one-step" 



The trace is over Dirac indices only and the subscripts on 
the quark propagator, S\{k), and constituent masses, Mi 
and M-2, denote quark flavors. Quark flavor is also indi- 
cated by the subscripts q and Q, where a meson of type 
m has the quark flavor structure m = qQ and m m de- 
notes the meson mass. The corresponding isospin factor 
and quark-meson coupling constant are labeled by C™ 
and g m qQ, respectively, and are determined within the 
NJL model [JO] HTJ . The integrated elementary splitting 
function is obtained from the elementary TMD splitting 
function via integration over p±, that is, 

d™(z) = J d 2 P± d^(z,pl) = -fg 2 nqQ z 

x f d 2 p ± y\ + [(z - 1)M! + M 2 } 2 

X J (2^) 3 \p\ + z(z - l)Ml + zMl + (1 - z)rr4] 2 ' 

(9) 

The probability densities d™(z) are then obtained 
by multiplying a normalization factor so that 
S m Jo dzd q n {z) = 1. The isospin and momentum 
sum rules are then satisfied automatically [101 HTJ . 



1 We use the following LC convention for Lorentz 4-vectors 
(a+, a~ , a±), a = —m (a ± a 3 ) and a± = (a 1 , a 2 ). 



process in the NJL model, using light-cone coordinates^ 
The NJL model which we use includes only four point 
quark interactions in the Lagrangian, with up, down and 
strange quarks (see, for example, Refs. [4"3"rfi"5] for de- 
tailed reviews of the NJL model). In the present work, 
we use the notation introduced in our previous stud- 
ies gnun]. 



The elementary fragmentation function for quark, q, 
to emit a meson, m, carrying light-cone momentum frac- 
tion, z, and carrying transverse momentum, p±, is de- 
picted in Fig. [5j In the frame where the fragmenting 
quark has zero transverse momentum, but a nonzero 
transverse momentum component —p±/z with respect 
to the direction of the produced hadron [H [17] , the un- 
regularized elementary TMD fragmentation functions to 
pseudoscalar mesons are given by 



-p-/z) 2ir5((p-k) 2 -M 2 ) 

(8) 
I 

Previously, we employed the Lepage-Brodsky (LB) reg- 
ularization scheme to calculate loop integrals such as that 
in Eq. (|9| . This method puts a sharp cutoff on the invari- 
ant mass squared, M 2 2 , of the particles in the final state 
(see Refs. [T71 HOI S3 SS] for a detailed description as 
applied to the NJL-jet model). The maximum invariant 
mass of the two particles in the loop, A12, is determined 
by 

M* 2 < K\ 2 = (J A* + pl + y/Al + p%\ , (10) 

where p± and p 2 denote the masses of the particles in 
the loop and A3 denotes the 3-momentum cutoff, which 
is fixed in the usual way by reproducing the experimental 
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FIG. 5. Feynman diagram describing the elementary quark 
to hadron fragmentation functions. 
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FIG. 6. The normalized integrated splitting functions for a 
u quark to tv + and K + , calculated using LB and LB-DIP 
regularizations with the same light constituent quark mass of 
M = 0.4 GeV. 



pion decay constant. For a light constituent quark mass 
of M — 0.4 GeV, the corresponding 3-momentum cutoff 
is A 3 = 0.59 GeV. The strange constituent quark mass is 
determined by reproducing the experimental kaon mass, 
giving the value M s = 0.61 GeV and the correspond- 
ing quark-meson coupling constants are g vq Q = 4.23 and 

9KqQ = 4.51. 

In loop integrals containing two particles, we assign a 
light-cone momentum fraction x (of the initial particle's 
light-cone momentum) to the particle with mass /zi and 
consequently a light-cone momentum fraction 1 — x for 
the particle with mass mu 2 . Then, in the frame where 
the initial particle's transverse momentum is zero, the 
invariant mass of the two particles in the loop can be 
expressed as 



Mr 



vi+pj , (4 +P± 

x 1 — x 



(11) 



The relation in Eq. (10), when applied to the integral 
in Eq. (|9|, yields a sharp cutoff in the integral over the 
transverse momentum, namely 



In this work, we employ a slightly modified version of 
the LB regularization, which replaces the sharp cutoff 
of the invariant mass squared in the integrals, namely, 
by a dipole regulator: 



e(A? 2 



Mf 2 ) 



1 



[i + (m 2 2 /a 2 2 ) 2 ] 



(13) 



A physical motivation for this regularization scheme is 
that it gives a pion quark distribution that at large x 
behaves approximately as (1 — x) 2S for Q 2 = 16 GeV 2 , 
which is in good agreement with the recent reanalysis of 
Aicher et al. which finds (1 — a;) 2 34 at the same scale [47j . 
Using this dipole cutoff version of the LB regularization 
scheme (LB-DIP), we fix the model parameters by re- 
producing the experimentally measured hadronic prop- 
erties, such as /tt and the kaon mass to determine the 
cutoff as A 3 = 0.773 GeV and the strange constituent 
quark mass becomes M s = 0.59 GeV. The correspond- 
ing quark-meson coupling constants are g^ q Q — 4.24 and 
9KqQ — 4.52. The quark distribution functions calcu- 
lated with LB-DIP regularization satisfy both the num- 
ber and momentum sum rules and allow us to set the 
model scale at Q\ = 0.2 GeV 2 in the usual way by com- 
paring the evolved pion distribution function with that 
obtained from experiment. This procedure is discussed 
in detail in Ref. [ft)] . 

The plots in Fig. [6] clearly show that the z range of the 
splitting functions calculated using LB-DIP allows for a 
smooth continuation of the corresponding splitting func- 
tions calculated using LB regularization to the endpoints 
z = and z = 1. The plots in Fig. [7] present results 
for the fragmentation functions of a u quark to 7r + and 
K + using the LB-DIP regularization scheme. We use the 
QCD evolution code of Ref. [50] at next-to-leading order 
to evolve our model results from the scale Qo = 0-2 GeV 2 
to Q 2 = 4 GeV 2 . We find a slightly better description of 
the empirical parametrizations compared to our earlier 
work [301 HI] , especially in the region where z is close to 
1. Previously, artifacts of the LB regularization did not 
allow a good description in this domain. 



p\ < Vl = *(1 - z) 



-1 2 



-{\-z)ii\-z [ i\. (12) 



A consequence of LB regularization is that it re- 
stricts the corresponding regularized functions to a lim- 
ited range of z, namely < z m i n ^ z ^ z max < 1, where 
z m i n and z max are d etermined by imposing the condition 
V\ ^ in Eq. (12). These range limitations depend on 
the masses of hadrons and quarks involved. For example, 
the z limits are very close to the endpoints (z = and 
z = 1) for quark splitting functions to pions, but are fur- 
ther from these endpoints for heavier hadrons like kaons. 
The plots depicted in Fig. [6] show the limited range for 
the normalized splitting functions of a u quark to n + and 
AT + , calculated using LB regularization. 



IV. TMD QUARK DISTRIBUTIONS IN THE 
NUCLEON 

The TMD quark distributions in the nucleon are again 
determined by utilizing the NJL model. The nucleon 
bound state is described by a relativistic Faddeev equa- 
tion that includes both scalar and axial-vector diquark 
correlations, where the static approximation is used to 
truncate the quark exchange kernel [35]. The relevant 
terms of the NJL interaction Lagrangian are 



Ci = G s (V J 5 Ct 2 P a /) (V C-^tiPa if) 

+ G a (^^CnT 2 /3 A ^ T ) (> T C-^tzuPa ifr) , (14) 
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FIG. 8. Feynman diagrams which give the unpolarized TMD 
quark distribution functions in the nucleon. The single line 
represents the quark propagator and the double line the di- 
quark f-matrix. The shaded oval denotes the quark-diquark 
vertex function, obtained from a relativistic Faddeev equation 
and the operator insertion has the form 'y + 5(x—^zp) \ (1 ± T3). 
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FIG. 7. The integrated fragmentation functions for a u 
quark to tv + (upper) and K + (lower), calculated using the 
LB-DIP regularization with a light constituent quark mass 
of M — 0.4 GeV and evolved from the model scale to 
Q 2 = 4 GeV 2 . The results are compared to phenomenological 
parametrizations of experimental data from Ref. |48] (HKNS) 
and Ref. [15] (DSS). The shaded area represents the uncer- 
tainties in the HKNS results. 



where C = 17270 and Pa — y I^a (j4 £ 2,5,7) are 
the color 3 matrices [35] ■ The strength of the scalar and 
axial-vector diquark correlations in the nucleon are de- 
termined by the couplings G s and G a , respectively. To 
regularize the NJL model for the calculation of the nu- 
cleon, we choose the proper-time scheme, with an in- 
frared and ultraviolet cutoff, labeled by Am and Ajyy, 
respectively. This scheme enables the removal of un- 
physical thresholds for nucleon decay into quarks, and 
hence simulates an important aspect of confinement |51l - 
153] . This simulation of quark confinement has also been 
shown to provide a natural saturation mechanism for nu- 
clear matter in the NJL model 53J. 

The proper-time regularization scheme is not used for 
the fragmentation functions because the emitted hadrons 
are not confined. Therefore, the confining nature of 
the proper-time regularization is not appropriate in this 
case. However, for consistency between both regulariza- 
tion schemes we use the same light constituent quark 
mass and fix the UV cutoff so as to reproduce the pion 
decay constant. 



The five parameters of our NJL model for the nucleon 
are the light constituent quark mass, M, the regular- 
ization parameters Am and Auv-i and the couplings G s 
and G a - These are determined by fixing M = 0.4 GeV, 
Am = 0.24 GeV, and then reproducing the nucleon mass, 
pion decay constant, and the nucleon axial coupling via 
Bjorken sum [37] [54]. Strange quarks are not yet included 
in our model for the nucleon. 

The leading-twist spin-independent TMD distribution 
of the quarks of flavor q in the nucleon is defined via the 
correlator l6l l55l 



Q(x, for) = P + J 



(2tt) 3 



x (N,S\^ q (0) 7 +W(OMr^T)\N, S) 



t+=o 



(15) 



where W(£) is a gauge link connecting the two quark 
fields, which are labeled by ip q . In QCD this gauge link 
is nontrivial for £x =^ 0, however at the level of approx- 
imation that we are working at, this gauge link equals 
unity in the NJL model. Our states are normalized using 
the noncovariant light-cone normalization, namely 



]T {N,S\$ q (0) 1 + i; q (0)\N,S)=3. 



(16) 



Equation ([15J) can be expressed in terms of two TMD 
quark distribution functions, namely, 



Q(x,kr) = q{x,k%) 



p — hi Q3 
£ T T _L / 7,2 

yr 1lT\ X ^ 



T )i 



(17) 



where the first TMD PDF integrated over fcx gives the 
familiar unpolarized quark distribution function and the 
second TMD PDF, known as the Sivers function [56ll57j . 
is time-reversal odd and is zero at the level of approxi- 
mation included in this work. 

To determine the TMD quark distributions in this 
model, it is convenient to express them in the form 



q(x, k T ) = —i 



dk+dk- 



k + 



Tr[ 7 +Af 9 (p,fc)] 



(18) 
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FIG. 9. Results for the u (upper) and d (lower) TMD quark 
distributions in the proton. 



where M q (p, k) is the quark two-point function in the 
bound nucleon. Therefore, within any model that de- 
scribes the nucleon as a bound state of quarks, the quark 
distribution functions can be associated with a straight- 
forward Feynman diagram calculation. 

The Feynman diagrams considered here are given in 
Fig. [HJ where the first diagram represents the so-called 
quark diagram and the second the diquark diagram. The 
single line in each diagram represents a quark propagator 
which is the solution to the gap equation and the double 
line is the diquark i-matrix obtained from the Bethe- 
Salpeter equation. The vertex functions represent the 
solution to the nucleon Faddeev equation. The resulting 
distributions have no support for negative x and therefore 
this is essentially a valence quark picture. By separating 
the isospin factors, the spin-independent u and d TMD 




FIG. 10. The Bjorken x dependence of (fey). Diquark corre- 
lations in the nucleon give rise to the quark flavor dependence. 



quark distributions in the proton can be expressed as 



u(x, k T ) — fq/x{x, k T ) + I fq/N{ x i fyr) 



+ 1 f 



q(D)/N\ 



^k 2 T ) + ir q(D)/N {x,e T ), (19) 



d(x, k T ) — | fq/ N (x, kj~) 

+ 5 fq{D)/N( x > fcr) + g fq(D)/N( x > for)- (20) 



The superscripts s and a refer to the scalar and axial- 
vector terms, respectively, the subscript q/N implies a 
quark diagram and q{D)/N a diquark diagram. Explicit 
expressions for the functions in Eqs. ( [T9| and (20 1 are 
given in the Appendix. 

Results for the u- and d-quark TMD quark distribu- 
tions functions in the proton are illustrated in Fig. [9] 
The Q 2 scale to which these results correspond is not 
determined by the model. Previously for the familiar 
spin-independent PDFs we fitted the valence u-quark 
distribution in the proton to the empirical result at 
some large Q 2 scale, this gives a model scale of Qq — 
0.16 GeV 2 [33 EZl EH] m the proper-time regularization 
scheme. Rigorous comparison with the experimental 
data requires QCD evolution of the model TMD PDFs, 
which is left for future work. Here, we just show the 
results as they emerge from our model, the exact scale 
of which is not so important for this purpose. When 
QCD evolution is included, both the TMD PDF and 
TMD fragmentation function model scales must be equal 
when determining observables like SIDIS cross-sections. 
The integral of these TMD PDF results over fcy gives 
the familiar spin-independent quark distributions func- 
tions, which satisfy the baryon number and momentum 
sum rules. The Bjorken x and k 2 -, dependence in these 
expressions is not separable, and therefore the Gaussian 
ansatz for the TMD quark distributions, namely, that 
they can be written in the form 



q(x, kr) = q(x) 



(21) 



7T (k 2 ) ' 

is not possible for our TMD PDF results. The Bjorken 
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FIG. 11. Results for the TMD M-quark distribution in the 
proton for fixed x slices, where the upper plot has x = 0.4 
and the lower plot is for x = 0.8. Also plotted are individual 
fits to the TMD quark distribution using the Gaussian ansatz 
of Eq.(21l for each x, with (fcy) the single fit parameter in 



each case. 



x dependence of (kj,) for our proton TMD quark distri- 
bution results is illustrated in Fig. [lOj where 



(k 2 T ) (x) 



J d 2 kx k\ q(x, k' 2 -,) 
J dPkx q{x, fcy) 



(22) 



If the x and k\ dependence of our TMD quark distribu- 
tions were separable then the curves in Fig. [10] would be 



constants, however we find that (fc 2 -) has about a 20% 
variation over the domain of Bjorken x. We also find that 
the x dependence of (fcy) for the u and d quarks differs 
somewhat, with the d quarks having slightly larger (kj.) 
for the majority of Bjorken x. 

Figure [TT] illustrates our TMD quark distribution re- 
sults at particular x values and compares them to a Gaus- 
sian ansatz fit for the same x slice. The Gaussian ansatz 
results are obtained by a least squares fit of the TMD 
factor in Eq. (21) to our ratios q(x, k 2 r)/q(x) calculated 
in the NJL model, using (k^) of Eq. (21 ) as the only fit 
parameter for each value of x. The fitted value of this 
parameter is approximately 20% smaller than the value 
of (fcy) calculated with our model distribution functions. 
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FIG. 12. TMD fragmentation functions for a u quark to tv + 
and 7r~. The upper figure illustrates the favored case, which 
peaks at relatively large z, while the unfavored case, shown 
in the lower figure, peaks at much smaller z. 



In the least squares fit, we in clud ed values of k 2 -, up to 
4 GeV 2 and the curves in Fig. 11 indicate that such a fit 
to a single Gaussian is reasonable only for a limited k\ 
region, for a single value of x. 



V. TMD FRAGMENTATION FUNCTION 
RESULTS 

In this section, we present NJL-jet model results for the 
TMD fragmentation functions. The number of emitted 
hadrons in the decay chain is set to Nn n k s = 6, which is 
sufficient to accurately obtain the pion and kaon fragmen- 
tation functions in the domain z > 0.02. We solve for the 
fragmentation of u, d, and s quarks to pions and kaons, 
utilizing Monte Carlo simulations and the expression in 
Eq. Q, similar to our previous calculations of the inte- 
grated fragmentation functions detailed in Ref . [41] . The 
computational challenge for the Monte Carlo simulations 
is to obtain sufficient statistics and this becomes signifi- 
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FIG. 13. TMD fragmentation functions for a u quark to K + 
and K~ . The upper figure illustrates the favored case, which 
peaks at relatively large z, while the unfavored case, shown 
in the lower figure, peaks at much smaller z. 



cantly more difficult when we include the transverse mo- 
mentum dependence, because now the number of bins be- 
comes quadratic in the size of the discrete bin size (taken 
to be 1/500 both for z and transverse momentum, in the 
corresponding units). Furthermore, the extent of the bins 
in the transverse momentum direction was extended to 
6 GeV 2 , in order to avoid any notable numerical artifacts 
arising from the limited range of transverse momentum. 
To overcome the numerical challenge, our software plat- 
form was developed to allow for parallel generation of the 
Monte Carlo quark decay cascades, with different seeds 
for their random number generators. The results were 
later combined to produce the high statistics solutions. 
The computations were facilitated on the small computer 
cluster at the Special Research Centre for the Subatomic 
Structure of Matter (CSSM) that consists of 11 machines 
with Intel Core i7 920 quad core CPUs running on the 
Linux Fedora Core 11 operating system and GCC 4.4. 
A typical calculation of fragmentation for a given quark 




FIG. 14. The averaged transverse momentum of ir and K 
mesons emitted by a u quark. 



type takes about 12 hours with 44 parallel processors. 

Results for the TMD favored and unfavored fragmen- 
tation functions for a u quark to it and K mesons are 
illustrated in Figs. [12] and [13] In each case, the favored 
TMD fragmentation functions have more support at large 
z, while the unfavored results are peaked at smaller z. It 
is also evident that the kaon fragmentation functions fall 
off more slowly in Pj_ than the corresponding pion frag- 
mentation functions. The drop in each of the fragmenta- 
tion functions for z < 0.02 is a consequence of choosing 
NLinks — 6, which means that in the Monte Carlo simu- 
lation there is a vanishingly small probability of emitting 
hadrons with z < 0.02. 

The Gaussian ansatz is widely used to describe the tra- 
verse momentum dependence of both quark distribution 
and fragmentation functions. In particular, the TMD 
fragmentation function of a quark q emitting a hadron h 
is often modeled by 



D h Jz,Pl) = D h Jz)- 



-PlKPl) 



(23) 



where D^(z) is the corresponding integrated fragmenta- 
tion function and (Pf) is the average transverse momen- 
tum of the produced hadron h, defined by 



Jd 2 P ± PlD^z,Pj) 
JcPP ± Dh(z,Pi) ■ 



(24) 



In analyses that assume a Gaussian ansatz for the TMD 
fragmentation functions, it is usual to assume that (P]_) 
does not depend on z, the type of hadron, h, or the quark 
flavor, q. These assumptions will be tested against the 
NJL-jet TMD fragmentation functions. 



The results in Fig. 14 depict the average transverse 
momenta of 7r and K mesons produced by a it-quark. 
These plots show that the average transverse momenta 
of the hadrons are relatively fiat versus z in the region 
0.3 < z < 0.6, however they have a significant depen- 
dence on the type of the hadron. We find that the av- 
erage transverse momentum of the kaons is significantly 
larger than that of the pions. 
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15. Normalized TMD fragmentation for the favored 
u — > 7r + process for z = 0.8 (upper) and the unfavored 
s — > K + process for z — 0.2 (lower). Also depicted are fits 
to the fragmentation functions using the Gaussian ansatz of 
Eq. (231, with (P±) as the single fitting parameter. 



The curves in Fig. [T5] depict the TMD fragmentation of 
a favored u — » tt + process for z = 0.8 and an unfavored 
,s — > K + process for z — 0.2. Also presented are least 
squares fits to the fragmentation functions for particu- 
lar z slices using the Gaussian ansatz of Eq. (23), with 
(P]_) the single fitting parameter for each z. The plots 
in Fig. [15] indicate that such a fit to a single Gaussian is 
reasonable only for a limited Pj_ region. Also, because 
(P]_) has a significant z dependence, the Gaussian ansatz 
for the entire TMD fragmentation function offers at best 
a crude approximation to the full results. The corre- 
sponding average transverse momenta obtained from the 
Gaussian fits are smaller than those obtained directly us- 
ing the relation in Eq. 
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VI. 



AVERAGE TRANSVERSE MOMENTA IN 
SIDIS 



For the SIDIS process we have created a Monte Carlo 
event generator that can calculate the physical cross- 
section. In future work this will enable us to analyze 
the relative importance of the different aspects of the 




FIG. 16. The averaged transverse momentum of ir + mesons 
in SIDIS produced on a u quark in a proton with light-cone 
momentum fraction x — 0.4. 




FIG. 17. The averaged transverse momentum of it and K 
mesons in SIDIS on a it quark in the proton with light-cone 
momentum fraction x = 0.4. The unfavored fragmentation 
functions rapidly approach zero for large z and this causes 
the dramatic changes in {P?) at large z illustrated in this 
figure. 



process and the implications of the constraints set in in- 
dividual experiments. We use it to determine the aver- 
age transverse momentum of the produced hadrons (at 
the model scale) observed in a SIDIS experiment, namely 
(Py), which is defined as 



/p2 w v Jd 2 P T P 2 D h q (x, Zl Pl) 
JdiP ± D%(x,z,p2) 



(25) 



The function DHx.z.P^) is defined in Eq. Q. The 
crosses in Fig. 16 represent results for (P^) acquired by 
7r + mesons in a SIDIS hadronization process, where the 
virtual photon strikes a valence u quark in a proton car- 
rying a light-cone momentum fraction of x = 0.4. We 
also plot as the dash-dotted line (Pf), which is the aver- 
age transverse momentum that the 7r + mesons acquire in 
the quark fragmentation process. Recall, that the trans- 
verse momentum P± is defined relative to the direction 
of the original fragmenting quark, while Pt is relative to 
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the direction of the photon momentum, these transverse 
momenta are related by Eq. ([IJ. For the factorization of 
the SIDIS cross-section given in Eq. ^ , it can be shown 
that (Pj.) is given by 



(P$)(x,z) 



(Pl)(z) 



z 2 (k^)(x) 



(26) 



As an additional check on the Monte Carlo calculation, 



in Fig. 16 we plot the result obtained from Eq. (26 1 as 



the solid line and find that it agrees perfectly with that 
obtained from the Monte Carlo event generator for the 
SIDIS cross-section. We also find that both (P]_) and 
(Pj) illustrated in Fig. 16 have a sizable z dependence. 



Illustrated in Fig. [17] are results for the average trans- 
verse momentum acquired by 7r and K mesons in the 
hadronization process in SIDIS, where the struck quark 
is a u quark in a proton with light-cone momentum frac- 
tion x — 0.4. The rapid approach to zero for the unfa- 
vored fragmentation functions in Fig. 17 is a consequence 



of the large z behavior of the unfavored (Pj) illustrated 
in Fig. [14J which also rapidly approach zero. The HER- 
MES experimental results for (Pj) measured in SIDIS 
on a deuterium target |27j , are of comparable size to our 
results shown in Fig. [it} We do not plot these HERMES 
results because the kinematic range is too different for 
a quantitative comparison. The average transverse mo- 
mentum of the kaons is larger than that of the pions at 
the low Q 2 scale of the model. Our model includes only 
the valence quarks in the proton, which should be the 
dominant component at x = 0.4. 



distribution functions [35], whose moments satisfy the 
baryon number and momentum sum rules. We found 
that the average fc T of the quarks in the nucleon have a 
significant x dependence and therefore the familiar Gaus- 
sian ansatz for the TMD PDFs produces only a crude 
approximation to our full TMD PDF results. 

Finally, using the TMD quark distribution functions 
for the nucleon and the results for the TMD fragmenta- 
tion functions, we constructed a Monte Carlo event gen- 
erator for the SIDIS process. Using this Monte Carlo 
event generator, we determined the average transverse 
momentum of the hadrons, (P^), produced in SIDIS. 
These results are of a similar magnitude to those ex- 
tracted from experiment, even at our relatively low model 
scale. As a cross check for this SIDIS Monte Carlo event 
generator, we compared our results for {Pj) with those 
obtained using Eq. (26), finding perfect agreement. We 



find that the (Pj) of the produced kaons is significantly 
larger than that of the pions, which is not apparent in 
the current experimental measurements. 

An interesting extension of our model would be to in- 
clude the vector meson and nucleon antinucleon emission 
channels. This extension has already been completed in 
our previous work on the integrated fragmentation func- 
tions. It would also be intriguing to consider the spin- 
dependent effects in the hadronization process, in par- 
ticular, to calculate the Collins fragmentation function. 
Further, using the NJL description of nucleon structure 
we will be able to develop a self-consistent description of 
the spin-dependent effects in SIDIS reactions. 



VII. CONCLUSIONS AND OUTLOOK 

In this work we extended the NJL-jet model to in- 
clude the transverse momentum dependence in the quark 
hadronization process. This was achieved using TMD el- 
ementary fragmentation functions and by keeping track 
of the quark's recoil transverse momentum in the hadron 
emission cascade. We modified the LB regularization 
scheme to remove artifacts that limit the z range of the 
splitting functions, and this in turn improved our de- 
scription of the integrated fragmentation functions. The 
TMD fragmentation functions for u, d, and s quarks to 
pions and kaons were determined using a Monte Carlo 
approach. The average P\ of the produced kaons was 
found to be significantly larger than that of the pions 
and in both cases (Pj_) had a sizable z dependence. The 
high statistical precision needed for these calculations 
was achieved through parallel computing on the small 
computer farm at CSSM. 

The TMD quark distribution functions in the proton 
were also determined using the NJL model. In this case, 
we used the proper-time regularization scheme, because 
this method simulates important aspects of confinement. 
Our TMD PDF results when integrated over fey give our 
earlier results for the familiar spin-independent quark 
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APPENDIX: NUCLEON TMD PDF 
EXPRESSIONS 

The u and d valence TMD quark distribution functions 
in the proton are given by 

u v (x, k T ) — fq/jj{x, k T ) + | fq/jsr{x, k T ) 

+ 2 fq(D)/N( x > kr) + I fq(D)/N( X ' &t)> (27) 

d v (x, k T ) = | fq/wix, kj) 

+ 2 fq(D)/N( x > kr) + g fq(D)/N( x i &y). (28) 

The individual quark diagrams terms have the form 
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/g/iv( X ) ^t) 
fq/N( X i ^t) 



16tt 3 
Z a Zn 
16tt 3 



(1 - x) I dr 
(l-x) 



1 + T X 



(M N + M) - m\ 



-T[k'^,+x(x-l)M%+xm 2 s + (l-x)M 2 



J dr [{al - 2a 2 a 3 - 2a%) (l+rx [(M N - M) 2 -ml])- 12a§ txMM n ] e~ T [ 
and the diquark diagrams terms are given by 



r \ki r +x{x~l)Mf l +xm 2 a + {l-x)M 2 



fq(D)/N(x, kr) = / / d v dz s ( x - v z ) fq/ D (z, M) f s q/N ( l - v) 



fq(D)/ N ( x i k r) = / / dydz 8{x - yz) f^ /D (z, k' T ) / g ° /JV (l - y). 



The diquark TMD quark distributions are 

3Z 



f S q /D^T) = ^ J dr [l + rx(l-x)m 2 s ] 



t \k^+x(x-l) m 2 + M 2 



47T 3 

f% /D (x,k 2 T ) = ^x(l-x) I dr [l+rx{l-x)ml 



2] e -T[k^+x(x-l)m 2 a +M 2 ] 



(29) 
(30) 

(31) 
(32) 

(33) 
(34) 



The constituent quark, scalar diquark, axial-vector di- 
quark, and nucleon masses in these expressions have the 
values M = 0.4 GeV, m s = 0.687 GeV, m a = 1.03 GeV, 
and Mjy = 0.94 GeV. The weight factors in the nu- 
cleon Faddeev vertex [351 GE3 EH] and its normalization 
are given by (ai, a 2 , a 3 ) = (0.429, 0.0244, -0.445) and 
Zn = 29.9, respectively. Finally the pole residues of the 
scalar and axial-vector diquark t-matrices are Z s = 14.5 
and Z a = 3.82, respectively. Integrating over kx in these 
expressions gives the familiar spin-independent quark dis- 
tribution functions, the moments of which satisfy the 
baryon number and momentum sum rules. 



These expressions have been derived using the proper- 
time regularization scheme, which in practice means to 
make the substitution 



1 



l 

(n - 1)! 



dr r n - l e 



rX 



(n- 1)! J 1/A 2 j 



1/A? fl 



(35) 



where X denotes the denominator function in a loop in- 
tegral after Feynman parameterization and Wick rota- 
tion. The infrared and ultraviolet cutoffs have the values 
A IR = 0.240 GeV and A uv = 0.645 GeV, respectively. 
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